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ABSTRACT 

In our previous paper we showed that the currently determined orbital pa- 
rameters placed four recently announced planetary systems HD 12661, HD 38529, 
HD 37124, and HD 160691 in very different situations from the point of view of 
dynamical stability. In the present paper, we deal with the last of these systems, 
whose orbital parameters of the outer planet are yet uncertain. We discover a 
stabilizing mechanism that could be the key to its existence. The paper is de- 
voted to the study of this mechanism by a global dynamics analysis in the orbital 
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parameter space related to the HD 160691 system. We obtained our results us- 
ing a new technique called MEGNO and verified them with the Fast Lyapunov 
Indicator technique (FLI). In order to be dynamically stable, the HD 160691 
planetary system has to satisfy the following conditions : (1) a 2:1 mean motion 
resonance, (2) combined with an apsidal secular resonance, (3) in a configuration 
P c (ap) — S — Pb(ap) (which means that the planets c and b may be considered 
as initially located at their apoastron around the central star S), (4) and spe- 
cific conditions on the respective sizes of the eccentricities. High eccentricity for 
the outer orbit (e c > 0.52) is the most probable necessary condition, while the 
eccentricity of the inner orbit e& becomes relatively unimportant when e c > 0.7. 
We also show that there is an upper limit for planetary masses (in the interval 
permitted by the undetermined line-of-sight inclination factor sinij) due to the 
dynamical stability mechanism. 

More generally, in this original orbital topology, where the resonance variables 
Q\ and #3 librate about 180° while 62 librates about 0°, the HD 160691 system 
and its mechanism have revealed aspects of the 2:1 orbital resonances that have 
not been observed nor analyzed before. The present topology with anti-aligned 
apsidal lines combined with the 2:1 resonance is indeed more wide-ranging than 
the particular case of the HD 160691 planetary system. It is a new theoretical 
possibility suitable for a stable regime despite relatively small semi-major axes 
with respect to the important masses in interactions. 

Subject headings: celestial mechanics, stellar dynamics - planetary systems - 
stars: individual (HD 160691) 

1. Introduction 

In our previous paper (Kiseleva-Eggleton et al. 2002) we applied the new technique 
invented by Cincotta & Simo (2000) and called MEGNO (the acronym of Mean Exponential 
Growth factor of Nearby Orbits) to a wide neighbourhood of orbital parameters determined 
using standard two-body Keplerian fits for the recently discovered multi-planetary systems 
HD 12661, HD 38529, HD 37124, and HD 160691 in order to distinguish between regular and 
chaotic planetary configurations. We showed that the currently announced orbital param- 
eters place these systems in very different situations from the point of view of dynamical 
stability. While HD 38529 and HD 37124 are located within large stability zones in the 
phase space around their determined orbits, the orbital parameters of the HD 12661 planets 
are located in a border region between stable and unstable dynamical regimes, so while its 
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currently determined orbital parameters produce stable orbits, a minor change within the 
margin of error of just one parameter may result in a chaotic dynamical system. The orbits 
in HD 160691 (Jones et al. 2002) at first appeared highly unstable, but using MEGNO we 
were able to identify a few stability zones in a parameter space which included the param- 
eters not determined from observations, such as the relative inclination i r between the two 
planetary orbits and the longitudes of the ascending nodes Vt. All these stable configurations 
are associated with the 2:1 mean motion resonance. The present paper is wholly devoted to 
a detailed and complete dynamical analysis of the HD 160691 planetary system by taking 
into account the angular orbital parameters not constrained by observational data (i r , Q), as 
well as sin i[ (the line-of-sight inclination factor) and the resulting different planetary masses. 
We have also explored the space of the mean anomalies M of the two planets taking initially 
into account the time of periastron passage r per given by Jones et al. 2002 (see Table 1). 

In one of our earlier papers (Gozdziewski et al. 2002) we have clearly identified (using 
MEGNO) the exact location of the 2:1 mean motion resonance and its width for the Gliese 
876 planetary system. A recent study by Hadjidemetriou (2002) of periodic orbits in this 
resonance predicts stable and unstable configurations of planetary systems depending on the 
hierarchy of planetary masses and eccentricities. The Gliese 876 system where m^, < m c 
and e h > e c is, according to this study, a stable configuration. In contrast, the hierarchies 
of the HD 160691 system are inverse, i.e. m b > m c and e& < e c . Hadjidemetriou (2002) 
found that a planetary system at the 2:1 resonance where the inner planet is much more 
massive than the outer planet is unstable for all values of the eccentricities. However, in 
the HD 160691 system the two planetary masses are comparable and this leaves the question 
about its stability open. 

In this work we have explored the parameter space available for planets in HD 160691 
in order to determine the stability conditions for this system. We have to notice that the 
orbital parameters of both planets are rather speculative due to the insufficient amount of 
observations, and even the existence of the second planetary companion is not yet fully 
confirmed (Jones et al. 2002, Butler et al. 2001). Nevertheless, the mechanism we present in 
this paper is probably the key to the existence of planetary systems like HD 160691. We have 
clearly identified the exact location of the 2:1 mean motion resonance, its width, and the 
secular resonance in apsidal longitudes preserving the stability related to the mean motion 
resonance. 
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Table 1. Orbital parameters of the HD 160691 planetary system (Data from Jones et al. 

2002, M* = 1.08 M Q ) 



Planet 


rripSmii (mj) 


a (AU) 


P (days) 


e 


lu (deg) 


Tper (HJD) 


b 


1.7 ± 0.2 


1.5 ± 0.1 


638 ± 10 


0.31 ± 0.08 


320 ± 30 


50698 ± 30 


c 


1.0 


2.3 1 


1300 


0.8 


99 


51613 



x In our numerical models we changed the value of a c from 2.3 AU to 2.381 AU (well 
within the error of its determination). The latter value gives the exact location of the 
2:1 resonance (by Kepler equation resolution) and is located in the middle of the stability 
valley on the [o^, a c ] parameter space (see Fig. la). 
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2. Method 

A classical method that allows us to distinguish between regular and chaotic dynamical 
states is the method of Lyapunov Characteristic Numbers (LCN). The estimation of LCN 
usually requires computations over long evolutionary times, sometimes much longer than the 
lifetime of the system studied. Let us note that chaotic in the Poincare sense means that the 
dynamical behavior is not quasi-periodic (the conventional definition useful for conservative 
dynamical systems), and does not necessarily mean that the system will disintegrate during 
any limited period of time. Let us state in addition that we use the property of stability in 
the Poisson sense : stability is related to the preservation of a certain neighbourhood relative 
to the initial position of the trajectory. In conservative systems, quasi-periodic orbits remain 
always confined within certain limits; in this sense they are stable. 

In the present work, we use two different methods in order to identify the dynamical 
state of the HD 160691 system : the MEGNO and FLI (the acronym of Fast Lyapunov 
Indicator) techniques. These two methods converge faster and are more sensitive than the 
LCN technique. 

FLI is the method introduced by Froeschle et al. (1997), permitting us to distinguish 
qualitatively between regular and chaotic motion in a dynamical system (see for example 
Pilat-Lohinger & Dvorak 2002). MEGNO is a new method developed by Cincotta & Simo 
(2000) that we have already successfully applied to the study of dynamical stability of ex- 
trasolar planetary systems (see e.g. Gozdziewski et al. 2001, 2002, Kiseleva-Eggleton et al. 
2002). This method provides relevant information about the global dynamics and the fine 
structure of the phase space, and it yields simultaneously a good estimate of the LCN with 
a comparatively small computational effort (Cincotta & Giordano 2000). The MEGNO is 
an alternative technique that proves to be efficient for investigation of both ordered and 
stochastic components of phase space (Cincotta, Giordano & Simo 2002). It provides a clear 
picture of resonance structures, location of stable and unstable periodic orbits, as well as a 
measure of hyperbolicity in chaotic domains (i.e. the rate of divergence of unstable orbits) 
which coincides with that given by the Lyapunov characteristic number. 
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Fig. 1. — Stability maps in the [a&, a c ] parameter space in 2-D (a) and 3-D (b) for the 
HD 160691 planetary system (resolution of the grid is 30 x 30 points). In (a), filled and 
open circles indicate stable orbits (< Y >= 2 ± 3%, and < Y >= 2 ± 5% respectively), 
while small dots (not surrounded by circles) indicate highly unstable orbits. < Y > is the 
MEGNO indicator characteristic value (Cincotta & Simo 2000). (a) is a cross section of (b) 
in the plane < Y >= 2 with projections of the different points according to their nature. In 
(b), the peaks indicate the magnitude of instability. The stability strip in (a) corresponds 
to a stability valley in 3-D (b) (however, its fine structure cannot be seen in this 3-D graph 
because of the scale factor in the z axis due to the high magnitude of the instability peaks). 
In each Figure of the paper, the parameters used (besides those scanned) have the following 
values (initial conditions for numerical integrations) : a b = 1.5, ej, = 0.31, % = 0°, Q b = 0°, 
oj b = 320°, M b = 156.3°; a c = 2.381, e c = 0.80, i c = 1°, fi c = 0°, u c = 99°, M c = 0°. 



Fig. 2. — Stability map in the [e&, e c ] parameter space for the HD 160691 planetary system. 
The symbols and the grid resolution are the same as in all Figures (see Fig. la). 



Fig. 3. — Stability maps in the [aj,, a c ] parameter space plotted in 2-D (a) and 3-D (b) as in 
Fig. 1 and computed for the HD 160691 system nominal elements (cf. Table 1 and Fig. 1) 
except for the mean anomaly of the b planet : M b = 0. Although the 2:1 resonance in mean 
motions is formally preserved, the stability valley is lost. 
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3. Stability conditions 

In order to identify the different dynamical behaviors in the parameter space, we use 
MEGNO which provides the exact location of stable and unstable orbits as well as a measure 
of hyperbolicity. Figures 1 and 2 show the dynamical state of the HD 160691 system as a 
function of both orbital semi-major axes a b and a c (Fig. 1) and eccentricities e b and e c (Fig. 
2). All other orbital parameters were taken from Table 1 with the addition of values of 
Q(b,c) and i(b, c ) undetermined from observations, and with r per replaced by the corresponding 
calculated values of mean anomalies M = (27r/P)(£ — r per ) (with an initial time of reference 
t = Tp er , one obtains M b = 156.3° and M c = 0°). In order to avoid a dynamical behavior 
bound to the plane, the two initial orbital inclinations % and i c are taken slightly different 
(0° and 1° respectively); due to gravitational interactions of the 3-body problem, the relative 
inclination i r is free to evolve in the 3-D space. In Fig. 2 as well as in other Figures in this 
paper a c was taken to be 2.381 AU according to the Table footnote. In all our Figures the 
intersection of horizontal and vertical lines indicates the 'observational' initial parameters 
taken from Table 1. By using the MEGNO indicator characteristic value < Y > (Cincotta 
& Simo 2000), filled and open circles in Fig. la (and all other Figures plotted in 2-D in this 
paper) indicate stable orbits (< Y >= 2 ± 3%, and < Y >= 2 ± 5% respectively), while 
small dots (not surrounded by circles) indicate highly unstable orbits (< Y > 3> 2). One 
can see that the system is globally highly unstable, except for a clearly marked stability 
region. Fig. la is a cross section of Fig. lb (plotted in 3-D) in the plane < Y >= 2 with 
projections of the different points according to their nature. In Fig. lb, the peaks indicate 
the magnitude of instability. The stability strip in Fig. la corresponds to a stability valley 
in 3-D (Fig. lb). However, its fine structure cannot be seen in this 3-D graph because of 
the scale factor in the z axis due to the high magnitude of the instability peaks. 

The stability strip is associated with the 2:1 mean motion resonance. Generally, the 
presence of a resonance structures the phase or parameter space in different regions, regular 
and irregular (i.e. chaotic). In Fig. lb, such different regions appear clearly. The regions 
on both sides of the stability valley present an irregular distribution of the values of < Y > 
meaning very different magnitudes of instability. The valley on the other hand presents a 
regular picture at the used scale, which is confirmed by the cross section (a) of (b). The strip 
indeed obtained in Fig. la, composed of filled and open circles, is sufficiently homogeneous 
and dense to be a stability zone. It is rather wide (~ 0.1 AU in both semi-major axes), and 
besides one order of magnitude larger than the one obtained for the Gliese 876 planetary 
system (also in 2:1 mean motion resonance, see Gozdziewski, Bois, & Maciejewski, 2002). 
However, for HD 160691 this stability valley may completely disappear for a small change in 
angular elements. The stability zone is indeed only permitted for a particular geometrical 
configuration of the orbits (determined by the combination of the elements Q, u, and i r ) 
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combined with particular relative positions of the two planets on their orbits determined by 
the M elements. For example, for a shift only on the element without losing the 2:1 mean 
motion resonance, close approaches at the periastron occur in the end, and consequently the 
dynamical behavior becomes totally and highly unstable as it is shown in Figure 3 (computed 
with Mb = 0° instead of 156.3°). The same situation appears with a comparable shift on the 
Cj elements. For instance, with ui> = 0, 100, or 200° instead of 320°, the stability valley is 
lost. The whole stabilizing mechanism depends on particular combinations of the elements 
f2, uj, i r , and M. We will discuss these stability conditions in detail below. 

Figure 2 shows two stability regions in the [e^, e c ] parameter space : (i) a small one with 
few points where both eccentricities are small, and (ii) a much larger region for high values of 
the outer orbit (e c > 0.5). Let us notice that in the latter, the eccentricity of the inner orbit 
eb is relatively unimportant when e c > 0.7. Let us emphasize that the large region with a 
sizeable number of stable points decreases significantly when a c is taken to be 2.3 instead of 
2.381 AU (see Fig. 3d in our previous paper Kiseleva-Eggleton et al. 2002). In the present 
Figure 2, due to the resolution of the grid, the dynamical behavior of the'observational' point 
in eccentricities (Table 1) does not appear. However, with a higher resolution (i.e. more 
integration points), computations confirm the existence of a stable point at the intersection 
of the horizontal and vertical lines. In addition, few unstable points are also very close. 
Let us recall that some zones (notably on the sides of stability structures) may present a 
sensitiveness to initial conditions for stable versus unstable points in a very close vicinity 
between them. This situation reminds us of the dynamical state of the HD 12661 planetary 
system (Kiseleva-Eggleton et al. 2002). 
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Fig. 4. — Stability maps in the parameter spaces [a>&, (a), i&] (b), [u c , i c ] (c), and 
[fib, h] (d) for the HD 160691 planetary system. The symbols and the grid resolution are the 
same as in all Figures (see Fig. la). 
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Figure 4 shows how the dynamical state of the HD 160691 system depends on the relative 
inclination i r between the planetary orbits and on the angular elements u and Q for both 
planets b and c. In the general case of a 2-planet system, the relative inclination i r is related 
to the longitudes of ascending nodes Qf, and Q c as follows : 

cos i r = cos if, cos i c + sin i b sin i c cos(f2 c — fib) 

However, in our computations, for a variation of % or i c , the inclination of the other orbit is 
taken equal to zero. As a consequence, i r = ib or i c , and then the inclinations on the vertical 
axes of the Figure 4 may be simply read as a scanning of the relative inclination i r between 
the orbits. It is important to remark that the distributions of structures are relatively similar 
in the pairs «&] and [u; c , i c ], as well as in [Q 6 , i 6 ] and [Q c ,i c ], with stable orbits appearing 
in strips of a comparable width ~ 100° around the nominal values of Ub and uj c (Table 1) for 
all values of i r . There are also additional stability zones, but very narrow, for small values 
of % r < 20° in both maps related to u, as well as in the two others related to Q. The nominal 
values for uj b and u> c , namely 320 and 99 degrees respectively, are located well inside the 
stable strips. 
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Fig. 5. — Stability map in the uj c ] parameter space for the HD 160691 planetary system 
(a) = Q + u). The symbols and the grid resolution are the same as in all Figures (see Fig. 
la). 



Fig. 6. — Stability map in the [M b , M c ] parameter space for the HD 160691 planetary system. 
The symbols and the grid resolution are the same as in all Figures (see Fig. la). 
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Figures 5 and 6 present the stability maps of the [ui,, Co c ] (where Cj = Q + uj) and 
[Mb, M c ] parameter spaces, respectively. Fig. 5 shows a linear relationship between the 
longitudes of periastron cDj, and uj c for stable configurations. This means that Ub and uj c 
precess at the same rate. Fig. 6 shows stable strips in the mean anomaly [Mj,, M c ] plane that 
permit us to identify adequate relative orbital positions for stable configurations at the 2:1 
resonance in mean motions. Figures 5 and 6 also show that the HD 160691 planetary system 
is dynamically stable only when the corresponding orbital parameters are simultaneously 
represented by stable points inside the strips in both [a>6,a) c ] and [M b ,M c ] planes. If only 
one of the two conditions is realized, the map [a&, a c ] exhibits a totally unstable dynamical 
state without any particular stability structure (see Fig. 3). We have also computed the 
[ub, oj c ] stability maps for different relative inclinations (from to 90 deg, not shown in this 
paper), and consequently we may point out that the whole mechanism of stabilization can 
be successfully applied even to high relative inclinations if Q,u,M are determined in an 
appropriate way : uj b and uo c on average precess at the same rate while the planets b and c, 
in 2:1 mean motion resonance, are located on their orbits in such a way that there are no 
close approaches at their periastron. The favourable respective positions are defined by the 
stable strips in the mean anomaly [M 6 , M c ] plane (cf. Fig. 6). 

The stability requires a particular geometrical configuration of the orbits (defined by 
elements Q, u>, and i r ), added to particular positions of the two planets on their orbits defined 
by the M element. In addition to the 2:1 mean motion resonance, the planets b and c have to 
be initially located in favourable respective positions, for example at their apoapses, in order 
to avoid close approaches at their periapses. Following analogous properties analyzed by Lee 
& Peale (2002) for the GJ876 planetary system, we have found in the HD 160691 system the 
simultaneous librations of the two mean motion resonance variables 6i = Xb — 2A C + cD& and 
6*2 = K — 2A C + Co c (where A = M + u), while the secular resonance variable 

$3 — $1 ~~ $2 = &b — &c 

librates about 180° with a period of about 2800 years and an amplitude of 80 degrees (see 
Fig. 7). 83 being only librating (i.e. without secular term), the two orbital planes on average 
precess at the same rate. This is a mechanism of secular resonance in periastron, which as a 
consequence ensures the maintenance of the orbital topology as well as the respective motions 
of the planets, that is to say in the present case, the 2:1 mean motion resonance without 
close approaches of the planets at their periapses. In other words, this mechanism, when 
combined with adequate relative positions of the planets on the initial osculating orbits, can 
be understood as the condition for the preservation of the dynamical stability related to the 
2:1 mean motion resonance. Because the two variables Q\ and 62 librate around 180° and 
0°, respectively (see Fig. 7), the lines of apsides are anti-aligned. This situation is different 
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from the GJ876 configuration in which the apsidal lines are aligned with 6±, 6 2 , and 8 3 all 
librating about 0° (Lee & Peale 2002, Gozdziewski, Bois & Maciejewski 2002). Let us note 
that in the familiar Io-Europa 2:1 resonance, the very small eccentricities lead to a geometry 
where 8\ is librating about 0°, while 82 and #3 are librating about 180° (Lee & Peale 2002). 
In this case, conjunctions occur when Io is near periapse and Europa near apoapse. 
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Fig. 7. — Librations of the mean motion resonance variables 9\ and 62, and of the secular 
resonance variable #3. The dynamical behavior of these three variables, derived from the 
integration of the geometrical elements given for the HD 160691 planetary system (see the 
values in caption of Fig. 1), express a mechanism of secular resonance in periastron with 
anti-aligned apsidal lines. 
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The whole stability mechanism for the HD 160691 system allows it to avoid close ap- 
proaches between planets, especially at their periapses (or near them). If it is not the 
case, the dynamical behavior is wholly unstable since the close approaches near periapses 
occur regularly (see Fig. 3 where Mb = M c = 0°). The original orbital topology of the 
HD 160691 planetary system associated to the 2:1 mean motion resonance can be written as 
PciflP) — S — Pb(ap) 1 , which means that the planets c and b may be considered as initially 
located at their apoastron around the central star S (state 1). Taking into account the 2:1 
mean motion resonance, the present topology is equivalent to S — P c (per) — Pb(ap) : after 
one revolution of the planet b, the planet c is at its periapse (state 2). This situation cor- 
responds (by taking into account the important amplitude of the secular resonance variable 
#3, namely 80 deg) to the nominal anomalies given in Table 1 (Mb ~ 160 deg, M c = 0). After 
one revolution more of the planet b, the planets are again near their apoastron (state 1). 
Besides, with the 2nd state, we understand that the present topology may drive to a solid 
stability where the sum of both eccentricities is particularly important (i.e. e& + e c > 0.7, 
high e c helps) (see Fig. 2). Because of the high eccentricities of the orbits, and despite 
relatively small semi-major axes, the relative distances between the two planets may remain 
sufficiently large over the whole evolutionary time scale of the system. 



1 This 3-D topology with anti- aligned apsidal lines is out of the scope of the various 2:1 resonance con- 
figurations analyzed by Hadjidemetriou (2002) (2 types of configuration) and Hadjidemetriou & Psychogiou 
(in preparation) (4 types in the planar case). 
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Fig. 8. — [at,, a c ] maps for different masses, k is the multiple of the determined minimum 
mass of the planets : (a) k — 1, (b) k = 1.3, (c) k = 1.4, (d) k = 1.5. The stable motion 
is represented by black squares (FLI < 10 8 ), "+" are plotted if 10 8 < FLI < 10 16 , and the 
white area shows the chaotic regions (FLI > 10 16 ). 
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We have also tested the robustness of this mechanism when the values of the planetary 
masses were progressively increased. In Figure 8, computed by the FLI technique in the 
[a,b, a c ] parameter space, the two masses are consecutively multiplied by a factor k — 1 (a), 
1.3 (b), 1.4 (c), 1.5 (d). The total destruction of the stable structure is reached for k = 1.5, 
which corresponds to a line-of-sight inclination factor sin 2; = 2/3 (i.e. i\ = 41.8°), as we 
know only the lower limit of planetary masses mo from the function mo = mp sin i\ where mp 
is the real mass. This result is in agreement with a numerical analysis by Kiseleva-Eggleton 
& Bois (2001) of multi-planetary systems of v And-type which showed that the dynamical 
regime, and consequently the lifetime of the system, depends strongly on the mass hierarchy 
as well as on the absolute values of planetary masses. In other words, assuming that the 
observed planetary system is dynamically stable, we could probably determine a window for 
the possible values of line-of-sight inclinations ii, and as a consequence give upper limits on 
the planetary masses. 



4. Conclusion 

Using the MEGNO technique of global dynamics analysis we scanned the most relevant 
cases of the orbital parameter space for the HD 160691 planetary system. We have found 
the existence of a stability zone ruled by a mechanism which involves angular elements of 
the system. This stability zone is indeed due to the 2:1 mean motion resonance coupled with 
adequate relative positions of the planets on their orbits avoiding close approaches at their 
periastron, the two apsidal lines being anti-aligned. The mechanism is not lost during the 
dynamical evolution of the system due to an apsidal secular resonance : the mean motion 
resonance variables are librating while the longitudes of periapse on average precess at the 
same rate. 

We conclude that in order to be dynamically stable, the HD 160691 planetary system 
has to satisfy the following conditions : (1) a 2:1 mean motion resonance, (2) combined with 
an apsidal secular resonance, (3) in a configuration P c (ap) — S — Pb(ap) (i.e. an apsidal 
antialignment), (4) and specific conditions on the respective sizes of the eccentricities. High 
eccentricity for the outer orbit (e c > 0.52) is the most probable necessary condition, while 
the eccentricity of the inner orbit e& becomes relatively unimportant when e c > 0.7. These 
four conditions, taking also into account various relative inclinations between the two orbits, 
determine the dynamical behavior of the system in such a way that the planets are never 
too close to each other. In the end, the HD 160691 system, where the resonance variables 
Q\ and 9 3 librate about 180° while 6 2 librates about 0°, has revealed resources of the 2:1 
orbital resonances that have not been observed nor analyzed before. The present orbital 
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topology [P c (ap) — S — Pb(ap)} combined with the 2:1 orbital resonance is indeed more wide- 
ranging than the particular case of the HD 160691 planetary system. It is a new theoretical 
possibility suitable for a stable regime despite relatively small semi-major axes with respect 
to the important masses in interactions. 2 

Combining our MEGNO maps (confirmed with the FLI method) for different pairs 
of parameters, it is possible to converge towards stability conditions related to different 
combinations of the angular parameters. As the observational determination of the elements 
of the HD 160691 system is far from being finalized, we hope that the maps presented here 
will be useful for testing both future observations and different parameter fitting techniques. 
In this respect, let us mention the fitting method of Laughlin & Chambers (2001), suitable 
for resonant interactions between the planets and where the true masses can be determined 
by eliminating the indeterminacy in sin %i inherent in fits that assume independent Keplerian 
motions. Let us note here that we have already tested some of the new fits for HD 160691 
obtained by Eugenio Rivera (private communication), and found that the only stable systems 
with k < 2 that do not meet our stability conditions would be ones with a very massive 
substellar distant companion (m c ~ 50m j) on a very large orbit (a c ~ 30 AU). This type of 
system is very different from the systems with Jupiter-mass planets on close orbits that we 
have discussed in this paper. 
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